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This manuscript discusses the computation of electromagnetic forces on a 
disc-type winding due to a standard switching impulse (SSJ). First, the 
resistances, inductances and capacitances (RLC) of a30 MVA, 33/11 kV disc- 
type distribution transformer were estimated to obtain the winding equivalent 
circuit. The transient voltage waveforms for each of the disc layers and 
corresponding resonances of the windings under the SSI were then obtained 
in time domains. Next, the axial and radial force distributions in the disc 
winding due to the SSI were computed. The forces on each disc layer and 
along the disc windings due to the SSI were computed based on the analytical 
and numerical methods via the finite element method (FEM) respectively. The 
non-uniform switching impulse voltage distribution results in non-uniform 
force distribution along the disc winding. The magnitude of the axially 


Transformer directed force on the disc winding is found to be higher as compared to the 
Transient radially directed force. 
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1. INTRODUCTION 

The overvoltages in the power system can be categorized into 2 types, which are external and internal 
overvoltages. The most common and severe external overvoltage is the lightning surge. The internal 
overvoltages generated from the electromagnetic interactions between its components are the switching and 
temporary surges. The presence of switching overvoltages in the transmission and distribution networks can 
affect the reliability of the transformers [1]. Switching overvoltage can be caused by routine operations, sudden 
changes in load, connection or disconnection of sections of the network and interruption of short-circuit current. 
It has long damped oscillating tail duration due to the electromagnetic transient and wave caused by the sudden 
changes in system operations [2]. The standard switching impulse (SSI) has a wavefront, Tp of 250 us and 
wave tail, T2 of 2500 us as per IEC 60060-1 [3]. 

The SSI can generate oscillations along the winding and causes abnormal electromagnetic stresses not 
only to the turn-to-turn but also in section-to-section winding insulation. It can also affect the linearity of the 
voltage distributions [4]. The magnitude of the switching impulses is approximately 4 to 5 times the nominal 
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system voltage [5]. Transformers should be able to withstand a series of electromagnetic forces during its 
operation. However, the cumulative stresses and forces due to the repeating switching surges on the winding 
can be severe and cause deviation in its electrical and mechanical characteristics. The electromagnetic forces 
are generated through the reactions between the leakage fluxes and current passing along the transformer 
windings. These forces are in axial and radial directions correspondingly to the leakage fluxes along the 
windings that are distributed in these components [6], [7]. 

The transients on transformer windings based on equivalent circuits and experimental models have 
been previously investigated in [8], [9]. Other studies have utilized the geometrical models of transformers via 
simulations and experimental setups to study similar effects [10]-[13]. These methods involve with a number 
of approximations due to the difficulties to obtain the transformer parameters since it involves with complex 
modelling of winding equivalent circuit. Multi-conductor transmission line (MTL) model approach was 
proposed to investigate the helical winding and continuous winding structures [14], [15]. The model involves 
the assumption that a pair of coils can be represented as a single capacitance. In recent years, electromagnetic 
models based on finite element method (FEM) have been introduced to support the modelling limitations on 
the geometrical and equivalent circuit models [16]-[19]. The transformer windings electromagnetic force 
measurements have been carried out in [20]—[23] with several simplifications and assumptions. FEM is found 
to be feasible for analysis of the electrical and mechanical damages on the windings [23]-[25]. FEM can model 
the actual capacitances between any two conductors in much more detail as compared to the MTL model. The 
FEM able to analyse the transient phenomenon based on the parameters that are computed from the transformer 
winding design information. In addition, the electromagnetic force of windings depends mainly on the 
transformer winding geometrical designs which in turn provides analysis closer to physical characteristics than 
the other modelling methods [26], [27]. 

Previous FEM studies mainly analyse the magnetic field distribution and electromagnetic force 
analysis to estimate the forces experienced by the windings under transient conditions. The winding model is 
quite general, which is only demonstrated as core-type transformers and constructed as layer windings [28]—[32]. 
General designs of low voltage (LV) and high voltage (HV) winding models were evaluated in [33]-[35]. 
Currently, there is limited work and understanding on the disc-type winding especially related to the 
electromagnetic force analysis to estimate the forces subjected to the winding. Comprehensive assessment is 
necessary to study the electromagnetic force distribution in the disc winding model design. This study can 
provide information on the electromagnetic force distributions on a winding model with detail disc winding 
configurations whereby the stress prone zone can be pinpointed. 

In this study, the electromagnetic force distribution on HV winding of a 30 MVA, 33/11 kV disc type 
transformer due to SSI is examined. The first section discusses the computation of the resistances, inductances 
and capacitances (RLC) elements of the HV disc winding, which is utilized as input data for the winding 
equivalent circuit model to obtain the transient resonances. The second section describes the simulation of the 
SSI waveform through a switching impulse generator circuit and it is applied as the excitation on the disc 
winding model to obtain the electromagnetic forces in FEM. The third section presents the computation of the 
electromagnetic forces on each disc layer based on analytical method and the overall electromagnetic forces 
subjected to the winding using the numerical method. The contribution of this study is the examination of 
electromagnetic forces phenomenon on a disc-type winding due to standard switching surge based on analysis 
of detail disc winding configurations. 


2. RESEARCH METHOD 

A 30 MVA, 33/11 kV Dyn11 transformer was examined in this study. The HV and LV winding’s 
nominal voltages are 34.5 and 13.8 kV. The HV and LV windings have BILs of 200 kV and 110 kV 
respectively. An assumption was made for the 33 kV winding whereby it was solidly grounded and the 13.8 
kV winding was low resistance grounded. The crest voltage of 165 kV for the SSI was set below the rated BIL 
of the HV winding with a protection margin of 21%. Only the 8 topmost HV disc windings were examined 
based on the assumption that the switching surge normally affected the outer HV winding. This study only 
discusses the HV disc winding since the magnetic flux that affects the electromagnetic forces and stresses 
primarily occurs on the corresponding winding. 


2.1. The standard switching impulse generation 

The impulse waveshape can be generated with a combination of a series RLC circuit under over damped 
condition. Figure | shows the SSI of 250/2500 us that was generated by using the impulse generator circuit [36]. 
The generator Simulink circuit is demonstrated in Figure 1(a) and the 165 kV magnitude SSI waveform is 
presented in Figure 1(b). The impulse generator circuit specification for the SSI is shown in Table 1. 
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Figure 1. Generating of SSI waveform through impulse generator circuit; (a) generator circuit and (b) 165 kV 
SSI voltage waveform 


Table 1. Impulse generator circuit specification for SSI 
Parameters Units 


The leakage inductance, L 0 
The primary capacitance, C; 10 uF 
The wave modulating capacitance, C2 1 pF 
The modulating inductance, Ln 0.9H 
The current limiting resistor, R 1000 Q 
The wave-tail resistor, Ry 1155Q 
The wave-head resistor, Ro 560 
400 V 


Amplitude of discharge voltage, Uci 


2.2. The HV winding RLC parameters calculations 
The single phase of the winding consists of 96 HV disc windings and 25 LV layered helical windings. 


The front cross-sectional view of a single-phase winding is illustrated in Figure 2. In addition, the geometrical 


designations for the HV winding are shown in Table 2. 
The RLC elements of the HV disc winding were computed through FEM via quasi-state Ansys Q3D 


and represented as RLC matrices according to [19]. The R matrix was established as the skin resistance of the 
winding conductor. The L matrix comprises of self and mutual inductances. The C matrix comprises of turn- 
to-turn, inter-disc and winding conductor to the ground capacitances. These parameters were computed from 
the structural specification of the winding using the ANSYS Q3D platform. The capacitances were established 
through method of moment (MoM) and fast multiple method (FMM) analyses via quasi-electrostatic solver. 
The AC resistance and inductance were estimated according to MoM via quasi-magnetostatic procedure. 
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Figure 2. The single-phase disc winding front cross-section view 


Table 2. The HV winding geometrical description 


Parameters Values 
Number of discs in one phase 96 
Number of turns per disc 30 
Height of conductor 11.5 mm 
Width of the conductor 2.4mm 
The thickness of the insulation (double-sided) 0.5 mm 
Distance between each disc 3mm 
Cooling duct between layers 12 and 13 5mm 
Inner radius of the HV winding 374.5 mm 
Outer radius of the HV winding 466.5 mm 
Total circumference of the HV winding 79.17 mm 
Height of HV winding 1437 mm 


Insulation between HV-LV windings (9 mm of oil + 1 mm of pressboard) 20 mm 


The turn-to-turn and inter-disc capacitances were computed based on MoM and FMM [19]. The turn- 
to-turn capacitances for every turn has been computed based on [37]. The mutual capacitance, Cm or inter-disc 
capacitance was calculated based on (1) [38]: 


27.6X1xey F (1) 
1012 in( [,-2#)’) 


where d is the distance between each disc and w is width of the conductor. €, is the relative permittivity of the 
insulation between the turns. The admittance matrix, [Y] = [G] + j@[C] which consisted of conductance [G] 
and capacitance [C] matrices was computed using the quasi-electrostatics solver [39]. The inductance (AC) 
and resistance (AC) were computed through MoM. The self-inductance, L,.,r was computed using (2): 
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al 


0.0021[In(—— 


)+0.5+0.2235(4*")| a 
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where /, h and w are the length, height and width of the winding conductor. The mutual inductance, Lm, exists 

within two winding conductors whereby one of the conductors is considered as the source that reacts with the 

neighbouring conductors upon a time-varying magnetic field via quasi magnetostatic interaction [38]. The 

impedance matrix was presented as Z = R + jwL. The magnitude of the mutual inductance was computed by 

(3) based on the coupling factor and conductors length [40]: 
_— 0.0021[1n(2+)-1+(4)] ia 


106 


(3) 


where / is the length of conductor and d is the space between the 2 conductors. Because of the skin and 
proximity effects, the resistance was represented as frequency-dependent and computed using (4) [41]: 


l 


R = ——Ohm (4) 


206(h+w) 


where / is the length of conductor, o is the conductors’ conductivity and 6 represents the conductor skin depth. 
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The impedance matrices were computed via Ansys Q3D parameter extractor with current as excitation using 
the quasi-magnetostatics setup. The computed RLC parameters result in 48x48 matrices. The RLC for each of 
the disc layer was extracted from the matrices as shown in Table 3. 


Table 3. RLC parameters of the 8 topmost HV disc winding 
Disc R (mQ) L (uH) Cu (pF) 
8  16.4483238 125.90481 8.12728 
19.5643928  62.032978 8.281119 
19.6474535 66.438095 8.280154 
19.6997534  322.64133 8.280349 
22.3123803  72.702426 8.28092 
15.4493239  39.693763 8.280631 
24.4277872 83.315436 8.280979 
19.6181349  66.35932 8.127188 


FPnwkuN WA 


2.3. The transient resonance due to the SSI 

Figure 3 shows the transformer RLC equivalent circuit to compute the transient resonances when 
subjected to the SSI. The transformer RLC equivalent circuit is presented in Figure 3(a). The schematic diagram 
of the 96 HV disc winding within the subsystem is shown in Figure 3(b). The analysis was carried out via 
transient resonance at the time of winding subjected to the switching impulse. The equivalent circuit model 
validation was performed by a comparison between the calculated and simulated voltage distributions [42]. 
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Figure 3. The transformer RLC equivalent circuit to compute transient resonances; (a) transformer RLC 
equivalent circuit and (b) subsystem of the 96 HV winding 
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2.4. The axial and radial forces based on analytical method 

The analytical method was employed to compute the forces subjected to each of disc layers. The radial 
F,, and axial F, directed forces were calculated and estimated according to the winding geometrical design 
according to (5) and (6) [43]. 


2x1? x(NI)2XDmX1077 


FE. = 7 


(5) 


FE, = 2x1? xAX(NI)*XDmX1077 (6) 
Neff 

where h is the length of the winding, herris the effective length of the radial flux path, calculated as 0.222 x h. 
Dy is the mean winding diameter, NI is the ampere turn of the winding and A is the winding's total length. The 
F,, was calculated based on the free distance between two windings as the axial flux density (Bz) between the 
windings. F,, was assumed to be constant along the length of the windings. It was generated by the integration 
of the instantaneous ampere-turn (NJ) in each of windings as well as the leakage flux derived as (5). The F, 
was calculated based on the residual ampere-turns method [44]. It is the relationship between the radial flux 
density (B,), the average of ampere-turns at the total length of winding and the effective length of radial flux 
path (he,--). The magnitudes of axial and radial forces in each disc layer with its leakage fluxes are presented 
in Table 4. 


Table 4. The axial and radial forces in each disc winding 
Disc Voltage dist.(kV) Fa(mN) Fr(mN) Ba(kT) Br (kT) 


8 164.89 88.029 30.858 0.103 1.857 
7 164.79 62.157 21.7889 0.087 1.560 
6 164.58 61.474 21.549 0.086 1.552 
5 164.37 60.990 21.380 0.086 1.546 
4 164.15 47.420 16.623 0.076 1.363 
3 163.94 38.163 13.378 0.068 1.223 
2 163.73 39.358 13.797 0.069 1.242 
1 163.52 60.863 21.335 0.086 1.544 


2.5. The electromagnetic force analysis based on numerical method 

In this case study, only 8 discs were modelled based on an actual geometrical design using Ansys 
Maxwell to analyse the electromagnetic force experienced by the winding due to switching transient 
phenomenon. Each of the disc consists of 6 conductors and 5 turns with high-frequency dependent 
computations. The HV winding geometrical details and RLC elements are shown in Table 2 and Table 3. The 
properties material of conductor, insulation and boundary were considered as copper, kraft paper and vacuum 
air. The complete model of the 8 topmost HV disc winding is displayed in Figure 4. 


Figure 4. 3D FEM model of 8-disc winding 


The electromagnetic solution was carried out in the transient analysis setup. The boundary condition 
was set as natural and the 165 kV SSI waveform was applied as the excitation on the winding to establish the 
transient phenomena. The electromagnetic force of the winding was then computed based on the Lorentz force 
analysis. The meshes were constructed by using adaptive length and skin depth-based refinement techniques 
to increase the accuracy of forces analysis in the winding. The winding model composes of 926,167 nodes and 
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369,168 elements. The transient was set from 1 us to 2500 us. In the 3D modelling, the calculation of the 
magnetic plane was defined as (x,y,z) since current was normal to the plane [45]. The component of leakage 
flux density was expressed as vector potential as in (7) [7], 


_ 10 (xAg) 
~ x Ox 


B, =-—* Baa: 8; (7) 
where B,, By and B_are directional components of leakage flux density and Ag is the magnetic vector potential. 
Lorentz force was employed to express the electromagnetic force on the winding as seen in (8) and (9) [6]: 


F=f) Jy x (Bk + B,A)dv = RL+EF (8) 
F, = B X Jy, Fy = B, X Jy (9) 


where, J,, is the y-directional current density, %, 7 and Z are the unit vectors in the winding plane. The F, and 
F, represent the radial and axial force, respectively. 


3. RESULTS AND DISCUSSION 
3.1. Transient voltage distribution under SSI 

The transient voltage distribution and resonances due to the SSI are shown in Figure 5. V96top is the 
applied impulse voltage while V96end, V93end, V83end, V73end, V63end and V53end are the voltages at the 
end of disc layers respectively. The transient resonances at the ends of disc layers indicate slight deviations 
between 0.41% and 9.04% at the tail time ranging from 350 us to 2500 us. The lowest and highest deviations 
are at V96end and V53end. As the number of disc layers increases, the deviations between V96top and the 
voltage at the end of disc layers increase, which indicates a decrement in the winding’s stress. The oscillations 
are influenced by the RLC components in the winding equivalent circuit. The transient voltage distributions 
along the disc layers initiate resonance oscillations, which causes the stress to build up in the winding. 
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Figure 5. Transient voltage distribution under SSI 


3.2. Axial and radial forces based on analytical method 

The axial and radial forces computed based on the analytical method via formulated calculations are 
presented in Figure 6. As seen in Figure 6, the magnitude of the axially directed force on the disc winding is 
greater than the radially directed force. This is anticipated due to the winding geometrical design. The disc 
winding structure causes a significant radial leakage flux component to emerge at the ends of the winding, 
which leads to the increment of the axial force. It is expected that the topmost disc layer, HV disc layer 8 
experiences the highest forces and resonances since it is the first layer that is subjected to the SSI. The forces 
acting on the other disc layers are expected to be uniformly distributed and lower in magnitudes as compared 
to the topmost disc. However, the forces of disc layers 4, 3, and 2 can decrease due to the RLC parameters in 
the ampere-turn element of the disc layers. 
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Figure 6. Axial and radial force in disc winding based on calculation 


3.3. Electromagnetic forces distribution based on numerical method 

The x, y, and z components of the electromagnetic forces based on the numerical method via electric 
transient analysis in FEM are presented in the time domain as seen in Figure 7. The axial force is seen 
significantly increases after 1.5 ms. The axial force is in positively directed force towards the HV winding, 
which represents a tensile force along the disc winding. The peak amplitude of the axially directed force of the 
HV winding is 35 mN. On the other hand, the radial force gradually increases to a negative amplitude after 1.5 
ms. The negatively directed force towards the HV winding indicates that it is an inwardly compressive force. 
The peak amplitude of the radially directed force on the HV winding is —2.2 mN. The HV disc winding 
experiences a higher magnitude of axial force as compared to the radial force. 

There are differences in the calculated and simulated forces due to the assumptions carried out in the 
analytical analysis and modelling stage for the numerical analysis. During analytical analysis, the thickness of 
winding insulation, mutual inductance inter-windings capacitance, spacers between windings, skin and 
proximity effect are not taken into account while these parameters are considered in the numerical analysis. 
The consideration of the spacers enhances the relative permittivity and thus increases the inter-disc 
capacitances for the numerical analysis. Nevertheless, the finding from both methods demonstrates a similar 
phenomenon which is the HV disc winding experienced higher axially directed force whether in each disc layer 
or in the overall disc winding when exposed to the SSI 
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Figure 7. The x (radial) and z (axial) force generated based on FEM 
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4. CONCLUSION 

The study shows that the non-uniform switching impulse voltage distribution results in non-uniform 
forces distribution along the windings. The axial forces experienced by each disc layer and on the HV disc 
winding are found to be higher in magnitude as compared to the radial force when subjected to the SSI. The 
winding electromagnetic force of disc winding is affected by the overall length of the winding as well as the 
configuration of the winding. This finding would assist the manufacturers on enhancing the radial winding 
structure during the design stage to minimize the axial force formation on disc-type winding. 
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